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Abstract 

The bifurcation diagram of a model nonlinear Langevin equation with delayed feedback is ob- 
tained numerically. We observe both direct and oscillatory bifurcations in different ranges of model 
parameters. Below threshold, the stationary distribution function p{x) is a delta function at the 
trivial state x = 0. Above threshold, p{x) ~ at small x, with a = — 1 at threshold, and 
monotonously increasing with the value of the control parameter above threshold. Unlike the case 
without delayed feedback, the bifurcation threshold is shifted by fluctuations by an amount that 
scales linearly with the noise intensity. With numerical information about time delayed correla- 
tions, we derive an analytic expression for p{x) which is in good agreement with the numerical 
results. 
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We obtain by numerical means the complete bifurcation diagram of a generic nonlinear 
and non Markovian Langevin equation that incorporates the effect of delayed feedback. 
Both pitchfork and Hopf bifurcation thresholds are observed and studied, and the results 
contrasted with two related limits: The deterministic limit of a differential delay equation, 
and the stochastic bifurcation of the same model without delay. 



The study of differential delay equations 



is an important topic in applied mathe- 



matics, with widespread applications in Physics (lasers, liquid crystals), control systems in 
Physiology (neural and cardiac tissue activity) 0, ^ , and Economy (agricultural commodity 
prices). Recent interest has arisen in the mathematical modeling of cellular function at the 
molecular level, especially in transcriptional gene regulation ^. Feedback regulation is a 
common motif in complex cellular networks, with delays arising from the complexity of the 
underlying network, or from the wide disparity in time scales of the many chemical processes 
involved in regulation For example, DNA is transcribed at a rate of 10 to 100 nucleotides 
per second, and it may take a delay of the order of minutes before the transcription factor 
appears as a finished product in the cell and is available for regulation. Significant delays can 
also be attributable to the time required for the diffusion of proteins through membranes, 
so that, for example, the auto regulated feedback on protein production at time t is often 
proportional to protein concentration at time t — r, where r is known as the delay time. For 
short delay times, a reaction may be approximated as being instantaneous, and the system 
as being in quasi equilibrium. However, when the delay is comparable to the characteristic 
time scales of reaction, the non instantaneous nature of the interactions can no longer be 
ignored, and delay terms need to be included in the governing equations for the network 
under study js, Q]. 

Experimental evidence has been mountirig that highlights the importance of stochastic 
effects in transcriptional regulation 
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, not only for natural networks, but for engi- 
neered gene circuits and networks as well [loi llll | . However, despite the wealth of evidence 
on the subject, delays in stochastic models of metabolic feedback are very often neglected, 
possibly because the resulting stochastic equations are no longer Markovian, and hence 
rarely tractable analytically. Exceptions include the derivation of a two time Fokker-Planck 
equation and the study of its small delay time limit in 12(], and results on the bifurcation of 
the first and second moments of a stochastic linear equation with delay j^, [l3]. We extend 
these latter results to the analysis of the stationary probability distribution function of a 



2 



nonlinear model, and discuss in detail the stability of the solutions that results from the 
interplay of delay and stochasticity. 

We focus on a canonical form of a nonlinear Langevin equation with multiplicative or 
parametric noise and delayed feedback 

X = ax{t) + bx{t - r) - x^{t) + ^{t)x{t) (1) 

where the constant a plays the role of the control parameter, b is the intensity of a feedback 
loop of delay r, and ^(t) is a white, Gaussian noise of intensity D. The initial condition is 
a function (pit) specified on t = [— r, 0]. We study the stationary probabihty distribution 
function p{x) for a range of values of a, b and D. 

The bifurcation diagram corresponding to Eq. ([1]) in the absence of noise is known (see, 
e.g. ^). Linearization around x = shows that trajectories decay asymptotically to zero if 

(I) b < -a and, (II) r < r, = ^7p= (2) 

The boundary separating exponentially decaying solutions from exponentially growing solu- 
tions is shown as the solid line in Fig. (jll). The upper branch, (I), is defined by Oc = —b and 
corresponds to a direct bifurcation (real eigenvalue), whereas the lower branch, (11), corre- 
sponds to a Hopf bifurcation (complex eigenvalue). In both cases, we show in the figure the 
case of r = 1, and hence the lower brach corresponds to Tc = 1 in Eq. ([2]). The cusp at the 
intersection of both boundaries is located at {a,b) = (1/r, — 1/r). 

The stochastic bifurcation analysis of Eq. ([!]) without delay (b = 0) is also known 
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|14l . 1151 . 116I | . Analysis of the linearized equation leads to the unphysical conclusion that the 
bifurcation threshold depends on the order of the statistical moment considered. With the 
saturating nonlinearity in Eq. ([1]), stationary probability distributions of x can be obtained 
both below and above threshold, and the location of the threshold properly determined. 
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The stationary distribution function of x with 6 = is 

a < —1 po{x) = 6{x) (3) 
a > -1 poix) = Nx'^e 2d (4) 

where the exponent a = a/D — 1, and is a normalization constant. The solution (jlj) exists 
but is not normalizable for a < — 1 and hence it is not a physically admissible solution. 
Therefore the bifurcation threshold is located at a = where p{x) changes from a delta 



function at the the origin to a power law at small x with an exponential cut off at large 
X. In this case, the existence of parametric fluctuations has no effect on the location of the 
bifurcation threshold: Both deterministic and stochastic equations exhibit a bifurcation at 
Oc = 0. In — 1 < a < 0, p{x) is unimodal with a divergence at x = 0, whereas for a > the 
distribution is bimodal reflecting saturation of x. 

We now turn to the case of delay, b ^ 0. Analytical results for the stability of the trivial 
solution X = of the linearization of Eq. ([1]) have been given in [s], jisl, and are shown in Fig. 
dl]). The bifurcation threshold of the first moment is shifted relative to the deterministic 
limit, only bounds have been given for the second moment and no results are available 
for p{x). Given the anomalous behavior described above for the stochastic bifurcation of the 
linear equation with 6 = 0, it is of interest to determine p{x) for the full model of Eq. ([T]). 
Unfortunately, the non Markovian character of this equation has precluded progress along 
these lines |12 |. 



We have first extended an existing high order algorithm for the integration of stochastic 
differential equations to the case of delay. The algorithm needs to take into account 
trajectories into the past for an interval r, and also new contributions from the stochastic 
terms that result from the coupling to the delayed feedback. In the numerical results shown 
below we employ, for technical reasons, an Ornstein-Uhlenbeck stochastic process ^(t) with 
intensity D and correlation time At, the same as the time step in the discretization of 
Eq. ([T]). The initial condition considered in all our calculations is a white and Gaussian 
random process in (— r, 0) of zero mean and intensity 1. The time step used is the numerical 
integration is At = 0.01. 

A qualitative view of the bifurcation is given in Fig. [H which shows a histogram of 
X once trajectories have reached a statistical steady state. For a ^ — 1, the histogram is 
approximately a delta function at x = 0. At a critical value Oc ~ — 1, the bifurcation point, a 
broad distribution emerges, although the most likely value remains x = 0. At larger values 
of a, the histogram becomes bimodal. This histogram shown corresponds to the direct 
bifurcation branch, but a similar graph is obtained for the Hopf bifurcation. Our results 
for the distribution function p{x) in these three ranges of values of a are shown in Fig. [2l 
For a < p{x) is approximately a power law with effective exponent a < — 1, but with a 
growing amplitude of p{x) at x ~ (not shown in the figure). Because of normalization, 
this growth implies a decaying amplitude for finite x, signaling a long transient leading to 
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the build up of the delta function at x = 0. Interestingly, the effective power law in the 
figure a ~ —1.2 < —1, indicating that would not be normalizable. For a > Oc we do 
obtain a time independent distribution with — 1 < a < 0. This distribution is normalizable, 
and represents the stationary distribution above threshold. We finally show j)[x) in the 
range of a for which the distribution is bimodal. The probability of the most likely value is 
approximately constant, but we still observe some transients in the vicinity of x = 0. Figure 
|3] shows the results of a power law fit to j){x) as a function of the control parameter a. We 
observe a smooth variation of the exponent a with a that allows a convenient determination 
of Oc, the value of a for which a = —1. This is the method that we have used to determine 
the bifurcation threshold in all the results presented below. 

We summarize our results for the bifurcation diagram of Eq. ([I]) in Fig. |H The analytic 
results for the threshold without noise = 0) are shown for reference, as well as the threshold 
of (x) of the linearized equation j^, and our numerical estimate. Except in the vicinity of 
the multi critical point (a, 6) = (1/r, — 1/r) our results are in excellent agreement with the 
analytic calculation of the linear equation, thus validating the accuracy of the numerical 
algorithm. As one approaches the point (1/r, —1/r) the Hopf frequency approaches zero 
and it is necessary to integrate the differential equation up to very long times to differentiate 
between an unstable trivial solution or an oscillation with a very long period. Since bounds 
on the threshold from (x^) for the linearized equation have been given in the literature |l3| . 
we also show in Fig. H] the threshold for this moment obtained numerically. Our numerical 
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results do agree with the known threshold for the special point of no delay 6 = given in |14j |. 
Interestingly, the thresholds for (x) and (x^) converge to the same values when one moves 
away from the point (1/r, —1/r), both in the direct and Hopf bifurcation branches. The 
figure also presents our results for bifurcation threshold defined directly from the stationary 
probability distribution function as discussed above. Our conclusion is that the stochastic 
threshold is shifted away from the deterministic threshold except in the special point of no 
delay (6 = 0), both for the direct and Hopf bifurcations. This threshold also agrees with 
that of the first and second moments of the linearized equation in the range of parameters 
in which both agree. Figure [5] shows the dependence of the shift in threshold as a function 
of the noise intensity D. In analogy with the case of no delay, we find a linear dependence 
in D. 

We next turn to the equation for p{x) that follows from Eq. ([T]). The difficulty in obtain- 
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ing a closed, analytic expression lies in the need to find the joint probability distribution 
p{x{t), x{t — t)). When r is larger than the correlation time of x, one can assume statistical 
independence between x(t) and xit — r), or p{x(t) , x(t — t)) =_p{x(t))p{x(t — t)) , an approx- 
imation that has been considered in the literature (e.g., ref. [6|). However, this assumption 
does not hold near a bifurcation since characteristic correlation times diverge. Instead we 
write p{x(t),x(t — r)) = p{x(t — T)\x(t))p{x(t)), where p{x(t — r)|x(t)) is the conditional 
probability of finding x{t — r) at t — r given x{t) at time t. We have then derived a the 
Fokker-Planck equation for p{x) that requires only the determination of {x(t — T)\x{t)), a 
correlation which we have not been able to compute analytically. We can, however, examine 
its behavior numerically. We have found that it reaches a stationary function (independent 
of t), and that for small x(t) is well described by — r)|x(t)) = (1 + a + 6 + D)x{t) (for 
r = 1). With this empirical relation and some straightforward algebra, we recover the same 
solution for p{x) given in Eq. (jlj) but with a = {a + h{l + a + h + D) — D) / D . This solution 
closely agrees with our numerical results of a versus a as shown by the solid line in Fig. [3l 
and with the value of the shift in Oc as a function of the noise intensity D shown in Fig. [5l 
In summary, many of the qualitative features of stochastic bifurcations under multiplica- 
tive noise are preserved under the addition of a delayed feedback. First, the bifurcation 
remains sharp. Since it is commonly assumed in the literature that the correlation time of 
X <^ r, the delay term in Eq. ([1])) would effectively act as an additive source of noise, 
leading perhaps to an imperfect bifurcation. We have shown this not to be the case. Second, 
and in agreement with the case of no delay (6 = 0), we observe that the moments of the lin- 
earized equation bifurcate at different values of the control parameter, which are themselves 
different from the threshold predicted from the distribution function of the full nonlinear 
equation. In contrast with the case of 6 = 0, however, the existence of delay introduces a 
shift in the bifurcation threshold, both for direct and Hopf bifurcations, shift that goes to 
zero as 6 ^ 0. The magnitude of the shift scales linearly with the noise intensity D. Finally, 
we have empirically derived the stationary solution of the distribution p{x) that agrees with 
our numerical determination of the bifurcation threshold. In view of our results, care must 
be exercised when analyzing bifurcation thresholds in numerical simulations of model gene 
regulatory networks when the analysis is based on the calculation of moments. 



6 



This research has been supported by NSERC Canada. 



[1] R. Driver, Ordinary and Delay Differential Equations, Vol. 20 of Applied Mathematical Sci- 
ences (Springer, New York, 1977). 
[2] M. Mackey and L. Glass, Science 197, 287 (1977). 

[3] L. Glass and M. Mackey, From Clocks to Chaos: The Rythms of Life (Princeton University 

Press, Princeton, NJ, 1988). 
[4] J. Hasty D. McMillen, F. Isaacs, and J. Collins, Nat. Rev. Gen. 2, 268 (2001). 
[5] M. C. Mackey and I. G. Nechaeva, Phys. Rev. E 52, 3366 (1995). 

[6] D. Bratsun, D. Volfson, L. S. S. Tsimring, and J. Hasty Proc. Natl. Acad. Sci. USA 102, 

14593 (2005). 

[7] M. Elowitz, A. Levine, E. Siggia, and P. Swain, Science 297, 1183 (2002). 
[8] N. Maheshri and E. K. O'Shea, Annu. Rev. Biophys. Biomol. Struct. 36, 413 (2007). 
[9] T. Kepler and T. Elston, Biophys. J. 81, 3116 (2001). 
[10] M. B. Elowitz and S. Leibler, Nature 403, 335 (2000). 

[11] J. Hasty, J. Pradines, M. Dolnik, and J. Collins, Proc. Natl. Acad. Sci. USA 97, 2075 (2000). 

[12] S. Guillouzic, I. L'Heurcux, and A. Longtin, Phys. Rev. E 59, 3970 (1999). 

[13] J. Lei and M. C. Mackey SIAM Journal of Applied Mathematics 67, 387 (2007). 

[14] R. Graham and A. Schenzle, Phys. Rev. A 25, 1731 (1982). 

[15] F. Drolet and J. Vifials, Phys. Rev. E 57, 5036 (1998). 

[16] M. San Miguel and R. Toral, in Instabilities and Nonequilibrium Structures, V, edited by E. 

Tirapegui and W. Zeller (Kluwer Academic, The Netherlands, 1999). 
[17] R. F. Fox, Phys. Rev. A 43, 2649 (1991). 



7 



-2 




-2-1012 



a 

FIG. 1: Long time histogram of x (in grey scale) as a function of the control parameter a from 
Eq. ([1]) with 6 = 1, r = 1, and D = 0.3. The histograms have been collected for t G (50,80) and 
further averaged over 60 independent runs. 
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FIG. 2: Probability distribution function p{x) for 6 = 1, r = 1 and D = 0.3 computed in the 
time interval t S (50, 80) and averaged over 10^ independent realizations. The values of the control 
parameter shown are: (top) a = —1.2, a ~ —1.25, (center) a = —1.1, a ~ —0.73, and (bottom) 
a = -0.9, a ~ 0.72. 
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FIG. 3: Fitted value of a as a function of a. We define the bifurcation threshold when a = — 1, 
or Uc — —1.15 for this parameter set (6 = 1, r = 1, and D = 0.3). The line is obtained from our 
empirical determination of the Fokker-Planck equation. There are no adjustable parameters. 
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FIG. 4: Bifurcation diagram of Eq. ([T]) with r = 1. The outer sohd hue corresponds to the 
deterministic hmit, the inner solid line is the analytic result of ref. [5] for (x), o and A, are the 
stability thresholds of (x) and (x^) of the linearized equation respectively, and * the threshold 
obtained from p{x). The three points labeled by • on the line 6 = are known results for no delay. 
All three are in good agreement with our numerical results. 



11 



CO 



-1.3 



-1.4 



-1 




FIG. 5: Bifurcation threshold Oc from p{x) as a function of noise intensity D. The parameters used 
are b = 1,t = 1, time averages for t G (300,350) and over 150,000 independent reahzations.The 
hne is obtained from our empirical determination of the Fokker-Planck equation. There are no 
adjustable parameters. 
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